Scaling behavior of optimally structured catalytic microfluidic reactors 
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In this study of catalytic microfluidic reactors we show that, when optimally structured, these 
reactors share underlying scaling properties. The scaling is predicted theoretically and verified 
numerically. Furthermore, we show how to increase the reaction rate significantly by distributing the 
active porous material within the reactor using a high-level implementation of topology optimization. 
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Chemical processes play a key role for the production 
and analysis of many substances and materials needed 
in industry and heath care. Generally, the optimiza- 
tion of these processes is an important goal, and with 
the introduction of microfluidic reactors involving lami- 
nar flows, the resulting concentration distributions mean 
better control and utilization of the reactors [|. These 
conditions make it possible to design reactors using the 
method of topology optimization Q, which recently has 
been applied to fluidic design of increasing complexity 

First, we report the finding of scaling properties of such 
optimal reactors. To illustrate the method we study a 
simple model of a chemical reactor, in which the desired 
product arises from a single first-order catalytic reaction 
due to a catalyst immobilized on a porous medium filling 
large regions of the reactor. 

Next, we show that topology optimization can be em- 
ployed to design optimal chemical micro-reactors. The 
goal of the optimization is to maximize the mean reaction 
rate of the micro-reactor by finding the optimal porosity 
distribution of the porous catalytic support. Despite the 
simplicity of the model, our work shows that topology 
optimization of the design of the porous support inside 
the reactor can increase the reaction rate significantly. 

Our model system is a first-order catalytic reaction, 

A — ► B, taking place inside a microfluidic reactor of 
length L, containing a porous medium of spatially vary- 
ing porosity 7(F) and a buffer fluid filling the pores. The 
porosity 7 is denned as the local volume fraction occu- 
pied by the buffer fluid [(| , and it can vary continuously 
from zero to unity, where 7 = is the limit of dense 
material (vanishingly small pores) and 7 = I is the limit 
of pure fluid (no porous material). The reactant A and 
the product B are dissolved with concentrations a and b, 
respectively, in the buffer fluid, which is driven through 
the reactor by a constant, externally applied pressure 
difference Ap between an inlet and outlet channel. The 
catalyst C is immobilized with concentration c on the 
porous support. 

The working principle of the reactor is quite simple. 
The buffer fluid carries the reactant A through the porous 
medium supporting the catalyst C. The reaction rate is 
high if at the same time the reactant A is supplied at a 



high rate and the amount of immobilized catalyst C is 
large. However, these two conditions are contradictory. 
For a given pressure drop Ap the supply rate of A is high 
if 7 is high allowing for a large flow rate of the buffer 
fluid. Conversely, the amount of catalyst C is high if 
7 is low corresponding to a dense porous support with 
a large active region. Consequently, an optimal design 
of the porous support must exist involving intermediate 
values of the porosity. Besides, the optimal design may 
involve an intricate distribution of porous support within 
the reactor, and to find this we employ the method of 
topology optimization in the implementation of Ref. |5j . 

In the steady-state limit, the reaction kinetics is given 
by the following advection-diffusion-reaction equation for 
the reactant concentration a, 



[11(7) ■ V] a = DV 2 a ~ k(y) a. 



(1) 



Here 11(7) is the velocity field of the buffer fluid, D is 
the diffusion constant of the reactant in the buffer, and 
— fc(7) a is the reaction term of the first order isothermal 
reaction, which depends on the concentration of the cata- 
lyst C through 7(r). In this problem three characteristic 
timescales t a , t r and t d naturally arise, 
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which correspond directly to the advection, diffusion, and 
reaction term in Eq. QJ, respectively. These time-scales 
will be used in the following theoretical analysis. Note 
that the index of f2 generally denote an average over the 
design region, e.g., kn = (Hl))n- 

The porosity field 7(1*) uniquely characterizes the re- 
actor design since it determines both the distribution of 
the catalyst and the flow of the buffer. In the Navier- 
Stokes equation, governing the flow of the buffer, the 
presence of the porous support can be modelled by a 
Darcy damping force density — a(j) u, where a is the 
local, porosity-dependent, inverse permeabilityQ. As- 
suming further that the buffer fluid is an incompressible 
liquid of density p and dynamic viscosity 77, the governing 
equations of the buffer in steady-state become 
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The coupling between a and 7 is given by the function 
a (7) = "max g g 1 +7 ) » where a max is determined by the 
non-dimensional Darcy number Da = - — ^-p- , and g is a 
positive parameter used to ensure global convergence of 
the topology optimization method|3,@- In this work Da 
is typically around 10 -5 , resulting in a strong damping of 
the buffer flow inside the porous support. The model is 
solved for a given 7(1") by first finding 11(7) from Eqs. i|3a|) 
and jSEjl and then a(r) from Eq. (J). 

Our aim is to optimize the average reaction-rate 
(A)(7) a)n of the reactor by finding the optimal porosity 
field 7(r). We therefore introduce the following objective 
function $(7), which by convention has to be minimized, 



$(7) = -(*(7)o)i 



(4) 



To better characterize the performance of the reactor 
and to introduce the related quantities, we first analyze 
a simple ID model defined on the x-axis. The porous 
medium is placed in the reaction region Q extending from 
x = to x = L. Eq. (|3bp leads to a constant flow velocity 
u, and as the complete pressure-drop occurs in the porous 
medium, we have p(0) = po + Ap and p{L) = pq. In this 
case the boundary conditions for the advection-diffusion- 
reaction equation Eq. Q arc a(— 00) = a , a'(— 00) = 0, 
and a' (00) = 0, where the primes indicate x-derivatives. 
We denote the outlet concentration a(oo) = ul- From 
Eqs. JTJ and Q, we then derive the following expression 
of the objective function 



$(7) =_(/c( 7 ) a } = - / [14(7)0' - £> a"] dx 

= <M(a L -a )-j[a'(L)-a'(0)]. (5) 

For simplicity, we now limit the analysis to the non- 
diffusive case (D = 0), and from Eq. (JSJ we get the objec- 
tive function defined in terns of the reaction conversion 
C, 

$ (7) = _I^ C; with C5 i_!£. (6 ) 
L ao 

With an explicit x dependence of the reaction rate co- 
efficient k(x), we obtain u(j)a' = —k{x) a with the solu- 
tion 

a(a-)=a e , with if(x) = / k(x)dx. (7) 

Jo 

This leads to the following expression of the conversion: 

C = l- e"^w =l-e"^w =l-e - ' 5 , (8) 

where we have introduced the dimensional-less 
Damkohler number 5 Hi 
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(9) 



having the physical interpretation of the ratio between 
the advection and the reaction timescale. 

To derive the flow speed u(^) in the ID model we first 
let q — > 00 resulting in 0(7) = a max (l — 7) and then by 
integrating Eq. I|3a|). = — p' — a max (l — 7) u, we get 
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(10) 



To solve the ID optimization problem analytically, we 
chose to abandon the spatial variations of 7 in the ID 
model. We have to find the solutions to 4^- = 0, and 
from Eqs. JHJ and (JSJ), we end up by having to solve the 
following equation 



l-e s + S(l + [3) =0, 



(11) 



where we have assumed that k(jn) cx (1 — jn) 13 ■ The spe- 
cific properties of the catalytic reaction determines the 
value of /?, e.g., if the full volume of the porous medium 
is active then (3 = 1, while if only the surface is active 
then (3 = 2/3. Solving Eq. (|ll(l gives the optimal value 
of 5f3, where the reference to [3 now is explicit. 

All numerical solutions are found using the commer- 
cial numerical modelling-tools MatlabQ and COMSOL 
[Tfij . To validate numerically the analytic results of the 
ID model, we solve Eqs. and i|l(J|) for a given homoge- 
neous design variable 70, and find th e op timal value using 
a brute- force optimization method |lll | . 

To obtain a general scaling parameter for the prob- 
lem defined in Eq. Q we reintroduce diffusion. How- 
ever, to minimize the trivial influence from the inlet and 
outlet, we only study the limit of low diffusion, e.g., 
t d t a ,t r . In this limit the optimal reactors involve 
a balance between the advection and reaction processes, 
and consequently we expect that t a and t r , should enter 
on equal footing in the scaling parameter. We are there- 
fore led to propose the following dimensional-less form of 
the scaling-parameter 
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Figure ^ shows that the measured values of 5p and 
Op for optimal porosity both scale with respect to 
\/ T A t r/ t d- The simulations cover 512 optimal reactors 
in a wide and dense parameter-scan [l2^ , and as they col- 
lapse almost perfectly on single curves, we have not dis- 
tinguished the data-points further. In the non-diffusive 
case D = and Jr A t r /t d = 0, exact values of 8p and Cp 
are determined by Eqs. JSJ and (| 1 1 p . and they match ex- 
actly with the numerical results, as seen in Fig.^ where 
they are marked by circles on the ordinate. 

We now introduce three types of 2D reactors: the 
uniform reactors, Fig. yfa), the membrane reactors, 
Fig. H[b), an d the topology optimized reactors, for which 
a few is shown in Fig. [21 First we optimize the simple re- 
actors in Fig. [21 They both depend only on one variable, 
which for the uniform reactor is the uniform porosity 7, 
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FIG. 1: Plot of ID results showing (a) Damkohler number 8 
and (b) conversion C, both as a function of ^/t a t r /t d and for 
optimal choices of porosity. In both cases = 2/3, 1 and the 
parameter-scan of each choice consist of 512 optimizations [l2T| 
which collapse nicely. For zero difffusion, ^/t a t r /t d = 0, 
Eqs. JIIJ and © give the exact results 82/3 = 0.9474, Si = 
1.2564, C 2 / 3 = 0.6122, and Ci = 0.7153, which are marked by 



circles on the ordinate axis. 



(a) Uniform reactor 



(b) Membrane reactor 
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FIG. 2: Illustration of the two simple 2D reactor setups, (a) 
the uniform reactor with porosity 7 and (b) the membrane 
reactor of width £ and with porosity 7* = 0. The horizontal 
dashed line is a symmetry-line. 
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FIG. 3: Representative collection of topology optimized reac- 
tor designs for deceasing values of yjT A t r /t d . (Left column) 
The distribution 7 of porous material in black together with a 
color-grading indication of the flow speed it. (Right column) 
The concentration a on top with the reaction rate k(-y)a be- 
low. Parameters (in SI units) L = 10~ 3 and the following 
values of [Da, D, Ap, k a ]\ (a) [10 -4 , 3 x 10" 8 , 0.25, 0.25], 



(b) [10~ 4 , 3 x 10~ 8 , 0.25, 1], (c) 
(d) [10 -4 , 10~ 8 , 0.25, 0.5]. 
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10~ 8 , 0.5, 1], and 



and for the membrane reactor is the width I of a porous 
membrane of porosity 7* = [Tll j. Because of mirror- 
symmetry in the xz-plane, only the upper half of the 
reactors are solved in all the following work. 

In the third type of 2D reactors we let the porosity 
7(r) vary freely within the same design region as for 
the uniform reactor. The optimal design is found us- 
ing the topology optimization method, described in de- 
tail in Ref. p|. This is an iterative method, in which, 
starting from an initial guess 70 of the design variable, 
the nth iteration consists of first solving the systems for 
the given design variable 7„, then evaluating the sensi- 
tivities by solving a corresponding adjoint problem, 
and finally obtaining an improved updated 7 n + i by use 
of the "method of moving asymptotes" (MMA) [ill Hi] . 
In Fig. |3 is shown a representative collection of topology 
optimized designs together with the corresponding flow 
speed it, concentration a, reaction rate k(j)a, and param- 



eter values. In the large parameter space under investiga- 
tion, our work shows a systematic decrease of pore-sizes 
and the emergence of finer structures in the topology op- 
timized reactors as the scaling parameter ^Jt a t r /t d is 
decreased. 

In Fig. 0] the conversion C is plotted as a function of 
\J T A t r/ t d f° r au optimal reactors of this work. It shows 
that all reactors collapse on curves similar to the ID reac- 
tors, although the topology optimized reactors exhibit a 
larger spread. We believe that this scaling is a signature 
of a general property of optimal immobilized catalytic re- 
actors. Note that the conversion of the uniform reactor 
in the low diffusion limit is a few percent higher than the 
theoretical estimate, an effect caused by low convection 
in the corners, resulting in 'dead zones'. 

In terms of the objective function $ the topologically 
optimized reactors are significant improved compared to 
the simple 2D reactors. To investigate the nature behind 
this improvement, we show in Fig. 0a log- log plot of the 
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FIG. 4: Overall scaling of the conversion C as a function of 
\/ t a t r/ t d f° r the different optimal reactors. The abscissa is 
logarithmic to emphasize the common scaling behavior. The 
dashed line indicate the theoretical value Ci for zero diffusion 
in the ID case. 
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FIG. 5: Log-log plot of the relation between convection C and 
flow rate Q for the different reactors, when normalized with 
values C UI1 if and Q un if of the uniform reactors. The reaction 
rate improvements are shown at the top (see text). 



flow rate Q and the conversion C normalized by the val- 
ues C un if and Qunif of the uniform reactors at the same 
parameters. Because Eq. (JHJ gives the following scaling 



of the objective function <& ~ QC, the rate of improve- 
ment with respect to the uniform reactors can be read 
off directly, as the contours of the improvement-factors 
of <& become straight lines, as showed by the dashed lines 
labelled by the corresponding factors in Fig. [5] It is seen 
that topology optimization can increase the reaction rate 
of the optimal reactors by nearly a factor 20, and further- 
more it does so by increasing the flow rate at the expense 
of lower conversions. The important insight thus gained 
is that the distribution of the advected reactant by the 
microfluidic channel network over a large area at mini- 
mal pressure-loss plays a significant role when optimizing 
microreactors. 

To conclude, we have analyzed a single first-order cat- 
alytic reaction in a laminar microfluidic reactor with opti- 
mal distribution of porous catalyst. The flow is pressure- 
driven and the flow through the porous medium is mod- 
elled using a simple Darcy damping force. Our goal has 
solely been to optimize the average reaction rate, with 
no constrains on the conversion or the catalytic prop- 
erties. A characterization of the optimal configuration 
has been derived theoretically and validated numerically. 
It shows an general scaling behavior, depending only on 
the reaction properties of the catalyst. The analysis is 
based on a very simple reaction since this emphasizes 
the points that the optimization of even simple reactions 
result in to non-trivial scaling properties and complex 
optimal designs. Using topology optimization to design 
optimal reactors give rise to reaction rate improvements 
of close to a factor 20, compared to an corresponding op- 
timal uniform reactor, and the improvement originates 
mainly due to an improved transport and distribution of 
the reactant. Furthermore, for the topology optimized 
reactors, we have found a systematic decrease of pore- 
sizes and the emergence of finer structures as the scaling 
parameter is decreased. Our work points out a new, gen- 
eral, and potentially very powerful method of improving 
microfluidic reactors. 
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